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The numerical implementation of an exchange-correlation functional is not always an accurate 
reproduction of its theoretical specification. For example, density functionals for exchange 
and correlation can be defined by an exchange or correlation hole function that is integrated 
with the Coulomb interaction to form an energy. This construction can be used to modify a 
density functional for use with any electron-electron interaction. Its most prominent use is in 
the Heyd-Scuseria-Ernzerhof (HSE) functional that generalizes the Perdew-Burke-Ernzerhof 
(PBE) model of exchange to a screened Coulomb interaction with an error function form. 
However, we find non-negligible numerical errors in the standard implementation of the HSE 
exchange hole integration. We formulate and implement a new method for evaluating the 
exchange hole integration that is simple, accurate, and efficient. Its numerical errors are 
bounded and minimized by applying basic elements of approximation theory. 



I. INTRODUCTION 

The Kohn-Sham density functional theory (DFT) of 
electrons is specified by density functionals that model 
the exchange and correlation energy of electrons. These 
functionals are presently limited by inaccuracies in ref- 
erence data, insufficient physical constraints, and sim- 
plification to convenient but approximate mathematical 
forms. In most cases, these functionals have an explicit 
analytic form and their numerical evaluation is straight- 
forward. However, increasing the accuracy of functionals 
inevitably requires more complicated forms. More com- 
plex functionals require improved numerical methods to 
evaluate them both accurately and efficiently. 

In this paper, we reassess the numerical evaluation of 
the semilocal exchange energy component of the Heyd- 
Scuseria-Ernzerhof (HSE) functional, which cannot be 
expressed exactly in a closed analytic form. Its evalua- 
tion requires additional numerical approximations and is 
subject to numerical errors. To date, it has had three dis- 
tinct implementations^"^. They differ primarily in their 
numerical integration of an electron exchange hole with 
a screened Coulomb interaction. They also embody dif- 
ferent regularizations of the underlying exchange hole 
model 4 , which breaks down in the limit of large elec- 
tron density gradients. We reevaluate these details and 
present a new method for approximating the integral 
that is simple and free of discontinuities and singular- 
ities. Our proposed reformulation is accurate to single 
precision (« 10 _T ) to reduce numerical errors far below 
the physical uncertainties of the HSE functional and limit 
the inevitable growth of errors in functional derivatives. 

It had been shown that the basic physics contained 
in the HSE exchange hole can be captured in more effi- 
ciently computable exchange hole models 5 . This is an ex- 
ample of using physical modeling to resolve the difficulty 



of numerically evaluating the HSE functional. However, 
those results are not numerically identical to the HSE 
functional and thus these constitute a distinct density 
functional. We consider an alternate approach: directly 
confront the numerical evaluation of the existing HSE 
functional as an application of approximation theory. 
This approach makes use of several basic elements of ap- 
proximation theory that include nonlinear minimax ap- 
proximation, special function evaluation, and error anal- 
ysis. Our result serves as a definitive numerical imple- 
mentation of the HSE semilocal exchange functional and 
demonstrates how increased sophistication in numerical 
analysis expands the set of numerically tractable models 
that might be used in the construction of future density 
functionals. 

The paper is organized as follows. First, we review 
the HSE functional. Second, we define a new regular- 
ization to ensure that its value and derivatives are well- 
defined. Third, we derive a new scheme for numerically 
integrating the exchange hole based on an approxima- 
tion of the complementary error function with Gaussians. 
Finally, we assess the accuracy of our new implemen- 
tation and compare with previous implementations for 
functional values and first derivatives. 



II. SUMMARY OF THE HSE FUNCTIONAL 

For clarity, we summarize the details of the HSE func- 
tional that are relevant to discussions of the exchange 
hole integration. The HSE correlation energy is identical 
to the Perdew-Burke-Ernzerhof (PBE) functional 6 and is 
not of interest here. The HSE exchange energy can be 
written as 
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The terms of this expression are defined by a splitting 
of the Coulomb interaction with error functions into a 
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short-range "screened" interaction and a long-range tail, 

1 erfc(W) erf(W) 

— = | . (2) 

The first term of Eq. ([!]) is the Fock exchange energy 
computed with the short-ranged screened Coulomb in- 
teraction, 
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which requires access to the spinful Kohn-Sham density 
matrix, /v,^' ( r ? r ')- The remaining terms replace the den- 
sity matrix with an exchange hole model 4 , J(r, |r — r'|), 
that is averaged over angle and assumes no mixing of 
electron spins, 

\p^(r,r')\ 2 « -2S^p a (r) 2 J(r,\r-r'\). (4) 

With an introduction of the specific PBE semilocal ex- 
change hole model, J PBE (s cr (r), kF^ a (r)\r — r'|), and a 
switch to inter-electron coordinates, u = r — r', the 
semilocal component of the exchange energy is written as 
the integral over a spatially local exchange energy density 
per electron, 

PCX) 

e^ E ' SR (r,w) =47r/v(r) / du uerfc(wu) (5) 
Jo 

xJ PBE (s a (v),k F , a (v)u) 
E™*^(u J )=Y: I drMr)e x P * E < SR (r,-). 

a J 

E^ be,lr (uj) has the same form, except with erfc(uu) 
replaced by erf (am). These terms depend on the elec- 
tron density, p a (r), the reduced density gradient, s a = 
| V p a \/(2kF i0 - Po), an d the Fermi wavevector, kp,a — 
(67r 2 pa) 1 ^ 3 • All equations and results in this paper are 
in Hartree atomic units. 

The J PBE function that defines the PBE exchange 
hole model is non-positive, obeys a normalization con- 
straint, 
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and reproduces the PBE exchange gradient enhancement 
factor, 



F PBE (s) = 1.804 
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1 + 0.27302857s 2 ' 
when integrated with the Coulomb interaction, 

dyyJ™ E (s,y) = -lF™ E (s). 
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With these constraints, J PBE was fitR based on properties 
of the uniform electron £ as, a principle of minimum in- 
formation, analytic integr ability of Eqs. (p| and d8j), and 



exact solubility of Eq. The exchange hole function 
only appears in calculations within integrals containing 
a Gaussian weight function, 
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with ci = 0.4516064, c 2 = 0.57786348, 

c 3 = 0.16520372, c 4 = 0.0068635965. 

Here g(x) contains the exponential integral function, 
Ei(x), 



g(x) := In 1 
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exp(x)Ei(x). (10) 



The variable / (denoted s 2 H(s) in the original papei^ 
is defined implicitly to satisfy Eq. (J8|, which now can 
be written as J(s, 0) = F PBE (s). The original version of 
Eq. ^ was in terms of a different set of five numerical 
coefficients, which we have rearranged into four different 
coefficients - Ci, c 2 , C3, C4 - for compactness and simplic- 
ity. 

The PBE exchange hole model used by the HSE func- 
tional has two mathematical pathologies that complicate 
its use. First, Eq. ([8| is satisfied by two values of / for 
8.26 < s < 11.14 and cannot be satisfied by any value of 
/ for s > 11.14. Second, all partial derivatives of Eq. ([9| 
with / diverge at / = (when 5 = 0) while all deriva- 
tives of F PBE (s) with s are finite. If f(s) is modeled as 
f(s) oc s m for s <C 1 (m = 4 was used in the original 
numerical fitpl) , the m th derivative of Eq. ^ with s will 
diverge erroneously at s = and cannot match the cor- 
rect response of F PBE (s). The exact f(s) is non-analytic 
at s = and approaches zero faster than a power law as 
s^0. 



III. REGULARIZATION OF THE EXCHANGE HOLE 

The primary goals of regularizing the exchange hole 
model are to satisfy Eq. ([8| for all s and to ensure all 
derivatives of the functional are well-behaved for all s. 
Secondary goals are simplicity and efficiency. 

We regularize the large-s limit without altering the 
small-s behavior by changing the explicit s-dependence 
in Eq. (l9|) from s to 
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which smoothly and monotonically limits the range of s 
from [0,oo) to [0,82). We make an arbitrary choice of 
Si := 8 and S2 := 11, which does not alter Eq. (JsJ) for 
small, physical values of s and guarantees that it has a so- 
lution for all values of s. This is similar to the previously 
proposed regularizatioiP, but strictly prevents unwanted 
alterations to the functional for s < s±. 

The regularized form of f(s) is defined to satisfy 
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The revised f(s) still causes numerical problems with 
derivatives of J(s, a) when the non-analytic behavior at 
s = is approximated by a rational function. This issue 
is resolved by either approximating f(s) with a form that 
exactly reproduces its non-analytic behavior at s = or 
by numerically evaluating f(s) as a solution of Eq. (12). 
We choose the latter option and apply Newton's method 
to Eq. (12) with f(s) initialized at its previous rational 
approximation^. 

All the singularities at / = in Eq. (JqJ) are con- 
tained in g(x) at x = 0. Its evaluation requires a careful 
cancellation of logarithmic singularities between the two 
terms. An algorithm^ for the numerical evaluation of 
Ei (x) prescribes a power series for small x and a contin- 
ued fraction for large x. We incorporate this result into 
two series expansions of g(x), for x < 2 as 
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where 7 « 0.5772156649015329 is Euler's constant, and 
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for x > 2. To converge g(x) to single precision for all x, 
each expansion requires ~ 20 terms. 



IV. NUMERICAL INTEGRATION OF THE EXCHANGE 
HOLE 

The difficulty of integrating the exchange hole with 
the range-separated Coulomb interaction originates from 
the error function in Eq. ([5|. The analytic form of the 
exchange hole was originally chosen by Ernzerhof and 
Perdew to simplify the integrals in Eqs. ([6| and ^ with- 
out considering the difficulty of other integrals. A simple 
method for integrating Eq. ([5| is to approximate it with 
something that is analytically integrable. Specifically, we 
seek to approximate the complementary error function as 
a sum of Gaussians, 



erfc(x) « Wj exp(— ajX 2 ). 
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Given such an approximation, we simplify the key inte- 
gral, 

- § / dy yerfc(by)J PBE (s, y)^T w ^ ( 16 ) 
Jo i=i 

Because the exchange hole is non-positive, we construct 
a simple bound on the error of this approximation, 
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with a relative error factor defined as 
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This analysis naturally limits the numerical errors in the 
HSE exchange energy to be a fraction of the PBE ex- 
change energy. Also, the error is independent of the scal- 
ing of the argument in erfc(x). 

Our numerical integration scheme is optimized by 
minimizing Eq. ( 18 ) over the choice of a$ and Wi, which is 
done once and then tabulated. This nonlinear minimax 
optimization is linearized to form 
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min max erfc(x) — 7 (w{ — ViX 2 ) exp(— a^x 2 ) 

Wi,Vi xe[0,oo) ^-^ 
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(19) 

This is solved using the Remez exchange algorithm, a 
standard tool in approximation theor}^. With these so- 
lutions, we continuously evolve the exponents as 5di = 
Vi/wi until convergence. This process is not globally 
convergent and requires a reasonable initial guess for a^, 
which is guided by solutions at smaller values of n. The 
goal of reducing the error factor e to below the single pre- 
cision machine-e, 1.2 x 10 _T , is achieved for n = 26 and 
the coefficients are given in Table [TJ There is no need to 
reduce e further because physical modeling errors in the 
HSE functional are orders of magnitude larger than the 
present numerical errors. 

With this numerical integration scheme, the HSE 
semilocal exchange energy density from Eq. ^ reduces 
to 
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TABLE I. Coefficients that minimize Eq. (18) for n = 26 
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accuracy brute-force numerical integration of the ex- 
change hole and the second HSE implementation 2 . While 
the second version of HSE is not the newest or most ac- 
curate, it is the most widely used. It is the only version 
available in vasfEEI espresscEED an d the present devel- 
opment version of LIBXC. The accuracy of the two im- 
plementations is shown in Fig. [I] We achieve the target 
of single precision relative error (« 10 -7 ) in the func- 
tional, while the old implementation has relative errors 
of up to 10" 2 . The errors are amplified by derivatives, 
but the stringent accuracy requirements of the new im- 
plementation keep first derivative errors below « 10 -6 . 
The relative error in the first derivatives gets as large as 
0.25 for the old implementation. 



VI. CONCLUSIONS 

An important long-term goal of electronic structure 
research is the continued improvement of accuracy within 
electronic structure simulations. In this paper, we have 
demonstrated an instance where numerical errors in a 
popular DFT exchange-correlation functional (HSE) are 
significant. We have analyzed these errors, devised an 
accurate and efficient method to eliminate them, and 
disseminated our new HSE implementation in the open- 
source LIBXC library. This result suggests that electronic 
structure research needs to increase its awareness of nu- 
merical errors and standards of numerical analysis. Im- 
provements in physical modeling of electron correlation 
alone cannot lower simulation errors below the floor set 
by numerical errors. 



V. ACCURACY VERIFICATION FOR FUNCTIONAL 
AND DERIVATIVES 



We have incorporated our new implementation of the 
HSE semilocal exchange functional into a development 
version of the LIBXC library of DFT functional^. This 
contribution includes first derivatives, which are mostly 
straightforward to derive. Derivatives of g(x) can be re- 
lated back to its value using the identity -jj^E\(x) = 



For example, its first derivative is 
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Derivatives of f(s) can be calculated from its implicit 
definition in Eq. (12), the first of which is 



l f ^ = i[ F x PBE W-W,o)] 

ds 



/(*) = 5s L ' x * r,'~, (22) 



We compare the new implementation to both a high- 
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FIG. 1. Accuracy of implementations of the HSE semilo- 
cal exchange energy. The new implementation is shown in 
the left column (a,c,e), and the old is shown in the right 
column (b,d,f). The first row (a,b) shows the energy den- 
sity per electron, £x,CT E ' SR ( r 5 The second row (c,d) shows 
the density derivative of the energy density per volume, 
^ y [p (T (r)e^ E ' SR (r,cj)]. The third row (e,f) shows the 
density-gradient derivative of the energy density per vol- 
ume, d | Vp t(r)| [^( r ) £ x,a E,SR (r,cj)]. The relative error be- 
tween high- accuracy values, X, and approximations, X, is 
plotted as . We set p to a typical material value of 0.01 

and consider a range of parameters: s £ [0,8] and uj G [0, 1]. 



